#################################
##' Code for:
##' Title: "Effective climate clubs require ambition, leverage, and insulation"
##' Author: Sam S. Rowan
##' Contact: sam.rowan [at] concordia.ca
##' Date: 2023-10-19
#################################
##' This script:
##' Insulation versus leverage in three sets of sectors
#################################


#### Initialize environment ####

## Working directory
# setwd("/Users/srowa/Dropbox/Projects/climate_clubs_trade")
# -- Set working directory to replication folder

## Requires objects from another script
source("code/06_ideal_type_clubs.R")
# comtrade_cbam (from 'd_comtrade_tidy.R')
# clubs_df (from 'e_ideal_type_clubs.R')
# wdi_emissions (from 'e_ideal_type_clubs.R')

## Load in trade data
comtrade_cbam <- read_csv("data/comtrade_cbam.csv")
comtrade_selected <- read_csv("data/comtrade_selected_carbon.csv")
comtrade_sampled <- read_csv("data/comtrade_sampled_sectors.csv")


#### Steel club ####

steel_eu_exposure <- comtrade_cbam %>% 
  # Keep only steel
  filter(bucket %in% "Steel") %>% 
  # Average steel flows over 2015:2019
  group_by(reporter_iso, partner_iso) %>% 
  summarize(trade_value_mean = sum(trade_value_mean, na.rm = T)) %>% 
  ungroup() %>% 
  # Get each country's world steel exports
  group_by(reporter_iso) %>% 
  mutate(exports_world = sum(trade_value_mean, na.rm = T)) %>% 
  ungroup() %>% 
  # Benchmark steel exports as share
  mutate(export_share_to_partner = trade_value_mean/exports_world) %>% 
  filter(partner_iso != "WLD") %>%
  select(reporter_iso, partner_iso, export_share_to_partner, exports_world) %>% 
  ungroup() %>% 
  as.data.frame() %>% 
  # Merge in EU membership data, twice: for reporter and partner
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("reporter_iso" = "iso3c")) %>% 
  rename(reporter_in_eu = eu_club) %>% 
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("partner_iso" = "iso3c")) %>%
  rename(partner_in_eu = eu_club) %>% 
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~as.numeric(.)) %>%
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~ifelse(is.na(.), 0, .)) %>% 
  # Now get trade flows
  rename(reporter_in_club = reporter_in_eu,
         partner_in_club = partner_in_eu) %>%
  # Create a dummy for whether an export relationship is into the club
  mutate(exporting_toclub = ifelse(partner_in_club==1,
                                   export_share_to_partner, 0)) %>%
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club) %>%
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>%
  ungroup() %>%
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, export_exposure_to_club, 
         export_value_world = exports_world) %>%
  distinct() %>%
  as.data.frame()

# Steel: scatterplot
p_steel1 <- steel_eu_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_in_club = as.factor(case_when(reporter_in_club==1 ~ "Member",
                                                reporter_in_club!=1 ~ "Non-member"))) %>% 
  mutate(export_world_log = log(1+export_value_world)) %>% 
  ggplot(., aes(x = export_world_log, y = export_exposure_to_club,
                label = reporter_iso, color = reporter_in_club)) +
  geom_smooth(method = "lm", se = F) +
  geom_text(size = 5) +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  scale_color_manual(values = c("#556B2F", "grey70")) +
  scale_x_continuous(breaks = c(0, 4.6, 9.21, 13.82, 18.4, 20.7, 23.02, 25.3),
                     labels = c("0", " ", "10K", "1M", "100M", "1B", "10B",  "100B")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(y = "Export exposure to club members",
       color = "Club member",
       x = expression(paste("Global steel exports, log transformed")))


# Steel: boxplot
p_steel2 <- steel_eu_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_in_club = as.factor(case_when(reporter_in_club==1 ~ "Member",
                                                reporter_in_club!=1 ~ "Non-member"))) %>% 
  ggplot(., aes(export_exposure_to_club,
                group = reporter_in_club,
                fill = reporter_in_club)) +
  geom_boxplot() +
  coord_flip() +
  scale_y_continuous(breaks = c(-0.2, 0.2), labels = c("Member", "Non-member")) +
  scale_x_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  scale_fill_manual(values = c("#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(x = "",
       fill = "Club member",
       y = "Club member") 

# Combo plot
cowplot::plot_grid(p_steel1, p_steel2, align = "v", ncol = 2, rel_widths = c(2/3, 1/3))
ggsave("text/comtrade_steel_eu_stacked.pdf", units = "in", height = 6, width = 12)

## Steel exposure to club members
steel_eu_exposure %>% 
  filter(reporter_in_club==0) %>% 
  arrange(-export_value_world) %>% 
  mutate(export_exposure_to_club = round(export_exposure_to_club*100, 2)) %>% 
  filter(reporter_iso != "CZE") %>% 
  # filter(reporter_iso != "GBR") %>% 
  slice(1:15) %>% 
  mutate(mean = mean(export_exposure_to_club))


#### Insulation and leverage for all trade, across sectors ####

## For whole economy, all sectors
total_eu_exposure <- comtrade_total %>% 
  # Merge in EU membership data, twice: for reporter and partner
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("reporter_iso" = "iso3c")) %>% 
  rename(reporter_in_eu = eu_club) %>% 
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("partner_iso" = "iso3c")) %>%
  rename(partner_in_eu = eu_club) %>% 
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~as.numeric(.)) %>%
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~ifelse(is.na(.), 0, .)) %>% 
  # Now get trade flows
  rename(reporter_in_club = reporter_in_eu,
         partner_in_club = partner_in_eu) %>%
  # Create a dummy for whether an export relationship is into the club
  mutate(exporting_toclub = ifelse(partner_in_club==1,
                                   export_share_to_partner, 0)) %>%
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club) %>%
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>%
  ungroup() %>%
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, export_exposure_to_club, exports_world) %>%
  distinct() %>%
  mutate(set = "Total trade") %>% 
  as.data.frame()

## Calculate leverage for all sectors
total_leverage <- total_eu_exposure %>% 
  filter(reporter_in_club==0) %>% 
  arrange(-exports_world) %>% 
  mutate(export_rank = row_number()) %>% 
  filter(export_rank <= 15) %>% 
  summarise(leverage_mean = mean(export_exposure_to_club),
            leverage_25 = quantile(export_exposure_to_club, 0.25),
            leverage_75 = quantile(export_exposure_to_club, 0.75),
            leverage_sd = sd(export_exposure_to_club)) %>% 
  mutate(bucket = "Total trade",
         set = "Total trade") %>% 
  relocate(bucket, .before = leverage_mean)

## Calculate insulation for all sectors
total_insulation <- total_eu_exposure %>% 
  filter(reporter_in_club==1) %>% 
  summarize(insulation_mean = mean(export_exposure_to_club),
            insulation_sd = sd(export_exposure_to_club),
            insulation_min = min(export_exposure_to_club, na.rm=T),
            insulation_10 = quantile(export_exposure_to_club, 0.10)) %>% 
  mutate(bucket = "Total trade",
         set = "Total trade") %>% 
  relocate(bucket, .before = insulation_mean)

#### Insulation and leverage for CBAM sectors ####

## Calculate exposure for CBAM sectors
cbam_eu_exposure <- comtrade_cbam %>% 
  # Merge in EU membership data, twice: for reporter and partner
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("reporter_iso" = "iso3c")) %>% 
  rename(reporter_in_eu = eu_club) %>% 
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("partner_iso" = "iso3c")) %>%
  rename(partner_in_eu = eu_club) %>% 
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~as.numeric(.)) %>%
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~ifelse(is.na(.), 0, .)) %>% 
  # Now get trade flows
  rename(reporter_in_club = reporter_in_eu,
         partner_in_club = partner_in_eu) %>%
  # Create a dummy for whether an export relationship is into the club
  mutate(exporting_toclub = ifelse(partner_in_club==1,
                                   export_share_to_partner, 0)) %>%
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club, bucket) %>%
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>%
  ungroup() %>%
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, bucket, export_exposure_to_club, exports_world) %>%
  distinct() %>%
  arrange(bucket, reporter_iso) %>% 
  mutate(set = "CBAM") %>% 
  as.data.frame()

## Leverage in CBAM
cbam_leverage <- cbam_eu_exposure %>% 
  filter(reporter_in_club==0) %>% 
  group_by(bucket) %>%
  arrange(-exports_world) %>% 
  mutate(export_rank = row_number()) %>% 
  filter(export_rank <= 15) %>% 
  summarise(leverage_mean = mean(export_exposure_to_club),
            leverage_25 = quantile(export_exposure_to_club, 0.25),
            leverage_75 = quantile(export_exposure_to_club, 0.75),
            leverage_sd = sd(export_exposure_to_club)) %>% 
  mutate(set = "CBAM")

## Insulation in CBAM
cbam_insulation <- cbam_eu_exposure %>% 
  filter(reporter_in_club==1) %>% 
  group_by(bucket) %>% 
  summarize(insulation_mean = mean(export_exposure_to_club),
            insulation_sd = sd(export_exposure_to_club),
            insulation_min = min(export_exposure_to_club, na.rm=T),
            insulation_10 = quantile(export_exposure_to_club, 0.10)) %>% 
  mutate(set = "CBAM")

  
## Plot total trade and CBAM
p_sectors <- cbind.data.frame(cbam_insulation, cbam_leverage[,2:5]) %>% 
  bind_rows(cbind(total_insulation, total_leverage[,2:5])) %>% 
  mutate(bucket = ifelse(bucket == "Total trade", "TOTAL", bucket)) %>% 
  mutate(set = ifelse(bucket == "TOTAL", T, F)) %>% 
  ggplot(., aes(x = insulation_mean, 
                y = leverage_mean, 
                color = set,
                fill = set,
                label = bucket)) +
  geom_point(aes(x = insulation_10, y = leverage_mean),
             pch = 19, 
             size = 2,
             color = "#556B2F") +
  geom_point(aes(x = insulation_min, y = leverage_mean),
             pch = 1,
             size = 1.5,
             color = "#556B2F") +
  geom_segment(aes(x = insulation_mean, 
                   xend = insulation_mean,
                   y = leverage_25,
                   yend = leverage_75),
               linetype = 1,
               color = "#556B2F") +
  geom_segment(aes(x= insulation_min,
                   xend = insulation_mean,
                   y = leverage_mean,
                   yend = leverage_mean),
               color = "#556B2F", 
               linetype = 2) +
  geom_segment(aes(x= insulation_10,
                   xend = insulation_mean,
                   y = leverage_mean,
                   yend = leverage_mean),
               size = 1.5,
               color = "#556B2F") +
  scale_fill_manual(values = c("white", "#556B2F")) + # "#E9F7EF"
  scale_color_manual(values = c("#556B2F",  "white")) + # "#E9F7EF
  geom_label(size = 6) +
  # geom_label(color = "#556B2F", size = 6) +
  labs(color = "Sector", x = "EU export insulation", y = "EU export leverage") +
  coord_cartesian(xlim = c(0,1),
                  ylim = c(0,1)) +
  theme_clubs +
  theme(legend.position="none")
ggsave(plot = p_sectors, filename = "text/comtrade_cbam.pdf", units = 'in', height = 6, width = 9)



#### High carbon, non-CBAM sectors #### 

## Comparing characteristics of CBAM sectors to non-CBAM sectors

## Selected high carbon sectors: exposure
selected_eu_exposure <- comtrade_selected %>% 
  # Merge in EU membership data, twice: for reporter and partner
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("reporter_iso" = "iso3c")) %>% 
  rename(reporter_in_eu = eu_club) %>% 
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("partner_iso" = "iso3c")) %>%
  rename(partner_in_eu = eu_club) %>% 
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~as.numeric(.)) %>%
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~ifelse(is.na(.), 0, .)) %>% 
  # Now get trade flows
  rename(reporter_in_club = reporter_in_eu,
         partner_in_club = partner_in_eu) %>%
  # Create a dummy for whether an export relationship is into the club
  mutate(exporting_toclub = ifelse(partner_in_club==1,
                                   export_share_to_partner, 0)) %>%
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club, bucket) %>%
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>%
  ungroup() %>%
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, bucket, export_exposure_to_club, exports_world) %>%
  distinct() %>%
  arrange(bucket, reporter_iso) %>% 
  mutate(set = "Other high carbon") %>% 
  as.data.frame()


## Selected high carbon sectors: leverage
selected_leverage <- selected_eu_exposure %>% 
  filter(reporter_in_club==0) %>% 
  group_by(bucket) %>%
  arrange(-exports_world) %>% 
  mutate(export_rank = row_number()) %>% 
  filter(export_rank <= 15) %>% 
  summarise(leverage_mean = mean(export_exposure_to_club),
            leverage_25 = quantile(export_exposure_to_club, 0.25),
            leverage_75 = quantile(export_exposure_to_club, 0.75),
            leverage_sd = sd(export_exposure_to_club)) %>% 
  mutate(set = "Other high carbon")

## Selected high carbon sectors: insulation
selected_insulation <- selected_eu_exposure %>% 
  filter(reporter_in_club==1) %>% 
  group_by(bucket) %>% 
  summarize(insulation_mean = mean(export_exposure_to_club),
            insulation_sd = sd(export_exposure_to_club),
            insulation_min = min(export_exposure_to_club, na.rm=T),
            insulation_10 = quantile(export_exposure_to_club, 0.10)) %>% 
  mutate(set = "Other high carbon")



#### Sampled, non-CBAM sectors #### 

## Getting exposure, leverage, and insulation for randomly sampled sectors

## Sampled sectors: exposure
sampled_eu_exposure <- comtrade_sampled %>% 
  # Merge in EU membership data, twice: for reporter and partner
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("reporter_iso" = "iso3c")) %>% 
  rename(reporter_in_eu = eu_club) %>% 
  full_join(., clubs_df %>% select(iso3c, eu_club), 
            by = c("partner_iso" = "iso3c")) %>%
  rename(partner_in_eu = eu_club) %>% 
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~as.numeric(.)) %>%
  mutate_at(vars(reporter_in_eu:partner_in_eu), ~ifelse(is.na(.), 0, .)) %>% 
  # Now get trade flows
  rename(reporter_in_club = reporter_in_eu,
         partner_in_club = partner_in_eu) %>%
  # Create a dummy for whether an export relationship is into the club
  mutate(exporting_toclub = ifelse(partner_in_club==1,
                                   export_share_to_partner, 0)) %>%
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club, bucket) %>%
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>%
  ungroup() %>%
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, bucket, export_exposure_to_club, exports_world) %>%
  distinct() %>%
  arrange(bucket, reporter_iso) %>% 
  mutate(set = "Sampled") %>% 
  as.data.frame()


## Sampled sectors: leverage
sampled_leverage <- sampled_eu_exposure %>% 
  filter(reporter_in_club==0) %>% 
  group_by(bucket) %>%
  arrange(-exports_world) %>% 
  mutate(export_rank = row_number()) %>% 
  filter(export_rank <= 15) %>% 
  summarise(leverage_mean = mean(export_exposure_to_club),
            leverage_25 = quantile(export_exposure_to_club, 0.25),
            leverage_75 = quantile(export_exposure_to_club, 0.75),
            leverage_sd = sd(export_exposure_to_club)) %>% 
  mutate(set = "Sampled")

## Sampled sectors: insulation
sampled_insulation <- sampled_eu_exposure %>% 
  filter(reporter_in_club==1) %>% 
  group_by(bucket) %>% 
  summarize(insulation_mean = mean(export_exposure_to_club),
            insulation_sd = sd(export_exposure_to_club),
            insulation_min = min(export_exposure_to_club, na.rm=T),
            insulation_10 = quantile(export_exposure_to_club, 0.10)) %>% 
  mutate(set = "Sampled") 

## Putting all of these summaries together
bind_rows(
  cbind(total_leverage[,1:2], total_insulation[,c(2,4,5)]),
  cbind.data.frame(bucket = "CBAM", 
                   leverage_mean = mean(cbam_leverage$leverage_mean),
                   insulation_mean = mean(cbam_insulation$insulation_mean),
                   insulation_min = mean(cbam_insulation$insulation_min),
                   insulation_10 = mean(cbam_insulation$insulation_10)),
  cbind(cbam_leverage[,1:2], cbam_insulation[,c(2,4,5)]),
  cbind.data.frame(bucket = "Other high carbon",
                   leverage_mean = mean(selected_leverage$leverage_mean),
                   insulation_mean = mean(selected_insulation$insulation_mean),
                   insulation_min = mean(selected_insulation$insulation_min),
                   insulation_10 = mean(selected_insulation$insulation_10)),
  cbind.data.frame(bucket = "Sampled",
                   leverage_mean = mean(sampled_leverage$leverage_mean),
                   insulation_mean = mean(sampled_insulation$insulation_mean),
                   insulation_min = mean(sampled_insulation$insulation_min),
                   insulation_10 = mean(sampled_insulation$insulation_10))) %>% 
  mutate_at(vars(leverage_mean:insulation_10), ~round(., 3)) %>% 
  select(bucket, leverage_mean, insulation_mean, insulation_10, insulation_min) %>% 
  knitr::kable(.,
               format = "latex", linesep = "",
               booktabs = T)
# -> paste into manuscript as table SI-1


sets_sectors <- bind_rows(bind_cols(cbam_insulation, cbam_leverage[,2:3]),
                          bind_cols(selected_insulation, selected_leverage[,2:3]),
                          bind_cols(sampled_insulation, sampled_leverage[,2:3])) %>% 
  as.data.frame() %>% 
  select(set, bucket, everything()) 


